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The process of relaxation of a system of particles interacting with long-range forces is relevant to 
many areas of Physics. For obvious reasons, in Stellar Dynamics much attention has been paid 
to the case of r~ 2 force law. However, recently the interest in alternative gravities emerged, and 
significant differences with respect to Newtonian gravity have been found in relaxation phenom- 
ena. Here we begin to explore this matter further, by using a numerical model of spherical shells 
interacting with an r~ a force law obeying the superposition principle. We find that the virial- 
ization and phase-mixing times depend on the exponent a, with small values of a corresponding 
to longer relaxation times, similarly to what happens when comparing for ./V-body simulations 
in classical gravity and in Modified Newtonian Dynamics. 
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1. Introduction 

Phase Mixing and Violent Relaxation are important dynamical phenomena, relevant in the astrophysical 
contest ( [Lynden-Bell, 1967], [Binney & Tremaine, 2008]). In fact ./V-body simulations revealed that the 
final states of gravitating systems experiencing a collapse from cold and clumpy initial conditions are 
structurally very similar to real elliptical galaxies ([van Albada, 1982], [Bertin et ai, 2005], [Nipoti et ai, 
2006]). In addition, the phase-space properties of the numerical end-products are remarkably interesting, 
with a differential energy distribution described very well by an exponential law, over a large value of 
energy values ([Binney, 1982], [Bertin & Stiavelli, 1984], [Ciotti, 1991]). Finally, recent investigations of 
Modified Newtonian Dynamics (MOND) showed that relaxation processes are much slower in MOND than 
in Newtonian gravity ([Ciotti et al., 2007], [Nipoti et ai, 2007]). This fact is relevant here, because in the 
weak limit MOND forces behave qualitatively as r -1 . Unfortunately, MOND is a non-linear theory, so that 
numerical simulations are more difficult than for systems with forces obeying the superposition principle. 
With this paper we begin an exploration aimed at understanding how the properties mentioned above 
depend on the specific nature of long-range (power-law) interactions. In particular, we study the relaxation 
of systems governed by additive r~ a forces, focusing on the simplified case of spherical systems. 
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2. The model 

Following previous studies ([Henon, 1962], [Takizawa & Inagaki, 1997], [Sanders, 2008], [Malekjiani et al, 
2009]), we numerically integrate the motion of N = 10 4 massive spherical shells. We assume that a surface 
element of each shell, of mass 5m and at vector position y, produces an acceleration at the point x given 
by 

5 = -G6m .. X ~, y j . 1 , (1) 

6 || x _ y|| a +! v ' 

where G is the "gravitational" constant, and the force index a, for reasons of convergence, is restricted 
to a < 3; the case of Newtonian gravity is obtained for a = 2, while a = — 1 corresponds to the case of 
harmonic force. The equations of motion for the shell i, of mass mi and radius r,, are 

dn dvi _ Jf A 

= Vi, — = di = -g + gii + ^ 9ji , (2) 

where Vi and Jj are the radial velocity and the (constant) angular momentum per unit mass of the shell, 
and ai is its effective radial acceleration. The field per unit mass acting on each surface element of the shell 
i due to the shell j at rj / n is 
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while 



9n = ~o -— s~ (4) 



is the self- field per unit mass acting on each surface element of the shell i; for a = 2 the celebrated Newton 
theorems are recovered. Note that ga can be obtained from gji by imposing nij = rrii and taking the limit 
Tj — > rf, the same limit is required to obtain the force in case of shell superposition. 

We performed experiments with different values of the exponent a, considering both radial collapses 
(Jj = 0) and collapses with Jj / 0. Numerically, when a shell collapses to the origin, the sign of the 
velocity Vi reverses to positive. As indicators of dynamical relaxation we consider the virial ratio 2K/\W\, 
the differential energy distribution n(E), and the phase-space evolution of the pairs (rj,^). We evolved 
the systems up to 50 dynamical times (idyn), where on dimensional grounds 



/ 2r a+1 

Mtot = J2i m i ls the total mass of the system, and ro is the half-mass radius at t = 0. We explored different 
radial distributions for the initial density profile, noticing that the final products are not strongly dependent 
on it. The initial conditions are characterized by a vanishing virial ratio, typical of cold collapses. The 
integration scheme ([Di Cintio, 2009]) uses a standard 2 th order leapfrog algorithm with constant timestep 
(e.g., see [Grubmiiller et al, 1991]). 



3. The virial ratio 

We first consider the time evolution of the virial ratio 2K/\W\. As the density distribution of the system 
can be written as 

N 

i=l i 



the virial function W R ( [Binney & Tremaine, 2008]) reduces to 

roc N ( 

W = 47r / p(r) g(r) r 3 dr = miri I 
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while the total kinetic energy of the system is 
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(7) 



(8) 



As apparent from Fig. 1, in purely radial collapses (and in collapses with Ji ^ 0, not shown), the evolution 
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Fig. 1. Top panel: evolution of the virial ratio over ~ 36fdyn f° r systems with — 0, and different values of the force index 
a (only bound shells have been considered in the plot). Low values of a correspond to larger and more long-lasting virial 
oscillations. Middle and bottom panel: kinetic energy and virial W in units of the initial binding energy U ln . Virial theorem 
and energy conservation show that for cold collapses in absence of escapers, W = 2(a — l)(7j n /(3 — a) at virialization. 



of the virial ratio shows a violent jump followed by an oscillatory phase, around the equilibrium value 1. 
The violent phase lasts ~ 2t(j yn , while the subsequent virial oscillations last considerably longer. A certain 
trend with a is apparent: as the exponent a increases towards 2 (the force range shortens), the peak of 
the first and of the successive jumps lowers. Remarkably, for the lowest value of a explored, (i.e., for the 
forces with the slowest decline with distance), the virial oscillations continue for several tens of dynamical 
times, similarly to what found in numerical simulations in MOND collapses. We interpret this result as a 
clear trend of the system behavior towards the harmonic oscillator case a = — 1, when from eqs. (3) and 
(4) it follows that cij = —GM tot ri + Jf/rf, so that each shell oscillates independently of the others and 
no energy exchanges between shells can take place. It is important to note that the somewhat arbitrary 
definition of tdyn in eq. (5) is found to be a quite accurate measure of the true dynamical time, as the first 
peak of the virial ratio occur at t ~ 2tdyn for all the explored values of a. 



For potentials that are homogenous functions of r, as here for a =fc 1, W = (a — l)U, where U is the potential energy. For 
q = 1 instead W = -GA/ t 2 ot /2 is constant, as for deep-MOND forces [Nipoti et al, 2007]. 
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4. Phase-space evolution 

The same behavior is also shown by phase-space evolution of the pairs (rj, Uj), represented in Fig. 2 for 
systems with a = 1.2 and a = 2. In the former case (top panels), the phase-space still presents an ordered 
structure after ~ lOidym symptom of a less efficient mixing than in the Newtonian case (bottom panels). 
This is reminiscent of the behavior of MOND systems, and the plots are strikingly similar to those in Fig. 2 
in [Ciotti et al, 2007]. This finding is also confirmed by experiments with a < 1 and a > 2 (not shown). 
From Fig. 2 (bottom panels) is apparent how shells are lost from the system with a = 2; note that ejection 
is energetically impossible when a < 1. 
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Fig. 2. From left to the right: phase-space at 0.7td yn , Ttdym am ^ ^*dya> f° r a system with a — 1.2 (top panels) and a = 2 
(bottom panels). The scale velocity is vt yp = rn/tdyni while rn and td yn are defined in Sect. 2. Note how phase mixing and 
shell ejection are more efficient in the case of Newtonian forces. 



5. Differential energy distribution 

We finally present the differential energy distribution at virialization, n(E). The quantity n{E)dE is defined 
as the number of shells with energy in the range (E, E + dE). The total energy of each shell is 

Ei = Ki + m % + <f>j^j . ( 9 ) 

where an integration gives the expressions for the potential cftji at rj due to the shell j, and the self-potentiaj^] 
<f>u. For a / 1 
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2 Note that —d4>ji/dri = while — d^j/dr, = 2(?jj. Note also that ^ is not the total energy of the system. 
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Fig. 3. Final differential energy distribution, normalized to E m - m /N, as a function of a. E m i D is the minimum energy value 
attained by the shell distribution. 



while for a = 1 

Gra ■ ( 1 \ 

<t>ji = ^ r ^[( r i + rj) 2 ln(ri +rj) - (fj - r.,) 2 In |rj - rj\ - 2^], (pa = Gnu ^lnr; + ln2 - -J . (11) 

The formulae for the self-potential can be obtained by direct integration or by taking rrij = nii and rj —> ri 
in the expression of <j)ji. The final n(E) for systems with 1.2 < a < 2 is shown in Fig. 3. The distributions 
(albeit quite noisy, due to the limited number of shells used in the simulations) roughly reproduce the 
expected decline with increasing (absolute) energy predicted by an exponential function 

n{E) = r M) e-^ E \, /3 > 0, (12) 

([Binney, 1982], [Ciotti, 1991]). However, the agreement of our results with eq. (12), known to provide an 
excellent fit over a large energy range of the virialized end-products of iV-body simulations of cold and 
clumpy systems, is quite poor. This is not surprising, as the imposed spherical symmetry of the shells does 
not allow the additional mixing effects due to non-radial density inhomogeneities. We conclude by noticing 
that the case a = 1.2 presents a curious, non-monothonic trend with energy; we verified that this feature is 
not a numerical artifact. We stress that also in MOND numerical simulations (that should correspond, in 
a broad sense, to the low-a cases) the final n{E) is very different to that of Newtonian simulations starting 
from the same initial conditions (see Fig. 5 in Nipoti et al, 2007). 



6. Discussion and conclusions 

With the aid of a simplified model of N massive, spherically symmetric and concentric shells, moving 
under the action of long-range, power-law additive forces, we studied the relaxation process leading to the 
final virialized states. We find that for all the explored values of a the relaxation process involves first a 
rapid Violent Relaxation phase, followed by a longer and gentle relaxation due to Phase Mixing. For a 
increasing from values < 1 to > 2, the relaxation time (in units of dynamical time) decreases steadily. In 
general, systems involving force laws with low a have longer mixing times compared to those of systems 
with high a, that are by contrast characterized by faster evaporation, resulting in a well marked core-halo 
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segregation. This finding can be useful to improve our understanding of MOND, where the relaxation 
mechanisms (due to the nonlinearity of the theory) are still quite poorly understood; in fact, the cases of 
low a present an evolution similar to that of iV-body MOND systems. Finally, the n(E) distribution of the 
final virialized states are qualitatively similar (with one exception), increasing for increasing energies, but 
are not described particulary well by an exponential law (especially for low a). The next step in our study 
will be the investigation of the relaxation of A-body systems with interparticle additive forces as in eq. (1). 
One of the main reasons for the successive analysis is to abandon the assumption of spherical symmetry, 
as it is expected that the increasing number of degrees of freedom will produce a faster relaxation, and 
a better agreement with an exponential n{E). Also, it will be interesting to investigate how radial orbit 
instability, characteristic of self-gravitating systems with most particles on radial orbits, depends on the 
specific nature of the force law. Unfortunately, for generic values of a, a field equation for the potential 
(such as Poisson equation for Newtonian gravity, or the p-Laplace equation for MOND systems) does not 
exists, so that particle-mesh methods cannot be used, and only the more time-expensive direct A^-body 
numerical schemes are available. An original code is currently under testing. 
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